##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--smg_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.molType.tsv",sep="")
    smg_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutsig_check/smg.list"
	images_path <- paste(work_dir,"/finalPlot/revise/smgs",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

mut_rate_gene_file <- opt$mut_rate_gene_file
smg_file <- opt$smg_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
    rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
    rgb(red=2,green=100,blue=190,alpha=255,max=255) 
    )


###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
smg <- data.frame(fread(smg_file))$Gene_Symbol
gene_order <- c("TP53","ARID1A","CDH1","APC","SMAD4","MUC6","PIK3CA",
              "CTNNB1","RHOA","ERBB2","CFTR","KRAS","MAP2K7","ARID2",
              "RNF43","TGFBR2","BMP6","FBXW7","CDKN2A","MTRR")

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol %in% smg)

dat_plot <- rbind( dat_mutRateGene )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IGC" , "DGC") , order = T )
dat_plot$value_text <- round(dat_plot$MutRate , 2) * 100 

sampleNum_record <- unique(data.frame(Class = dat_mutRateGene$Class , Molecular.subtype = dat_mutRateGene$Molecular.subtype , SampleNum = dat_mutRateGene$SampleNum))

###########################################################################################
## 计算P值
dat_plot$p.value = ""
dat_plot$p_text = ""

result_tmp <- c()

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

for(geneN in gene_order){
    print(geneN)
    for(molN in unique(dat_plot$Molecular.subtype)){
        tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Molecular.subtype %in% molN & Class == "IGC" )

        if(nrow(tmp_1)==0){
            tmp_1 <- tmp_2
            tmp_1$Molecular.subtype <- molN
            tmp_1$SampleNum <- as.numeric( subset(sampleNum_record , Molecular.subtype %in% molN & Class == "IGC")$SampleNum )
            tmp_1$MutNum <- 0
            tmp_1$MutRate <- 0
            tmp_1$value_text <- "0"
            tmp_1$Class <- "IGC"
        }

        tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Molecular.subtype %in% molN & Class == "DGC" )
        if(nrow(tmp_2)==0){
            tmp_2 <- tmp_1
            tmp_2$Molecular.subtype <- molN
            tmp_2$SampleNum <- as.numeric( subset(sampleNum_record , Molecular.subtype %in% molN & Class == "DGC")$SampleNum )
            tmp_2$MutNum <- 0
            tmp_2$MutRate <- 0
            tmp_2$value_text <- "0"
            tmp_2$Class <- "DGC"
        }

        tmp <- rbind( tmp_1 , tmp_2 )
        tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
        p <- fisher.test(tmp_fisher)$p.value

        if( p < 0.01 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
        }

        tmp$p.value <- p
        tmp$p_text <- ""
        tmp$p_text[1] <- p_text


        result_tmp <- rbind( result_tmp , tmp )
    }
}

result_tmp$Hugo_Symbol <- factor( result_tmp$Hugo_Symbol , 
    levels = gene_order , order = T)

###########################################################################################

result_use <- result_tmp

result_use$p <- result_use$p.value

#result_use$p_text=ifelse(result_use$p>=0.05,"","*")
#result_use$p_text=ifelse(result_use$p<0.05 & result_use$p>0.01,"*",result_use$p_text)
#result_use$p_text=ifelse(result_use$p<0.01 & result_use$p>0.001,"**",result_use$p_text)
#result_use$p_text=ifelse(result_use$p<0.001 ,"***",result_use$p_text)

result_use$p_pos <- 0.95
result_use$p_pos <- as.numeric(result_use$p_pos)

result_use$label_pos <- 0.5

result_use$percent_pos <- result_use$MutRate + 0.03
result_use$percent_pos <- as.numeric(result_use$percent_pos)

## 提取TP53
result_plot <- subset( result_use , Hugo_Symbol == "TP53" & Molecular.subtype %in% c("GS" , "CIN") )

#print(result_use)
result_plot$Molecular.subtype <- factor( result_plot$Molecular.subtype , levels = c("GS" , "CIN") , order = T )
result_plot$Class <- factor( result_plot$Class , levels = c("IGC" , "DGC") , order = T )
p <- ggplot(result_plot,mapping = aes(x = Class , y = MutRate , fill = Class)) +
    geom_bar(stat = 'identity', position = 'dodge' , width = 0.8 , color = 'black') + 
    facet_grid(.~Molecular.subtype , scales = "free")+
    theme_bw() +
    scale_y_continuous(
        breaks = seq(0,1,0.2) , 
        label = c( 0 , "20" , "40" , "60" , "80" , "100")
        ) +
    geom_text(aes(label= value_text , x = Class, y = percent_pos ), position=position_dodge(0.9),size=4, vjust=0 ,color="black", face='bold' , family="Helvetica")+
    geom_text(aes(label=p_text , y = 1.05 , x = 1.5),parse = TRUE,size=4.5 , color = "black") +
    #geom_text(aes(label=Type , x = 18 , y = label_pos ),size=5,family="Helvetica")+
    xlab("") +
    ylab('Percentage of samples with mutations (%)') +
    scale_fill_manual(values = col) +
    theme(
        title =element_text(size=4, face='bold'),
        legend.title = element_blank(),
        legend.text = element_text(size = 10),
        legend.key.width = unit(1, "cm"),
        legend.key.height = unit(1, "cm"),
        legend.position = 'none' ,
        strip.text =  element_text(size = 13 , color = 'black' ) ,
        axis.text.x = element_text(size = 10,color="black") ,
        axis.title.y =  element_text(size = 12 , color = 'black' ) ,
        axis.text.y =  element_text(size = 10 , color = 'black' , family="Helvetica") ,
        axis.ticks.length = unit(0.2, "cm") ,
        panel.grid=element_blank() 
    )

width <- 6
height <- 4
images_name <- paste0(images_path , "/MutRate.compare.TP53.GS_CIN.pdf")
ggsave( images_name , p , width = width , height = height )
images_name <- paste0(images_path , "/MutRate.compare.TP53.GS_CIN.tsv")
write.table( result_use , images_name , row.names = F , sep = "\t" , quote = F )
